arXiv:1505.00846v2 [cond-mat.quant-gas] 17 Sep 2015 


Dissipative fluid dynamics for the dilute Fermi gas at unitarity: 

Anisotropic fluid dynamics 

M. Bluhm and T. Schafer 

Department of Physics, North Carolina State University, Raleigh, NC 27695 

Abstract 

We consider the time evolution of a dilute atomic Fermi gas after release from a trapping po¬ 
tential. A common difficulty with using fluid dynamics to study the expansion of the gas is that 
the theory is not applicable in the dilute corona, and that a naive treatment of the entire cloud 
using fluid dynamics leads to unphysical results. We propose to remedy this problem by including 
certain non-hydrodynamic degrees of freedom, in particular anisotropic components of the pres¬ 
sure tensor, in the theoretical description. We show that, using this method, it is possible to 
describe the crossover from fluid dynamics to ballistic expansion locally. We illustrate the use of 
anisotropic fluid dynamics by studying the expansion of the dilute Fermi gas at unitarity using dif¬ 
ferent functional forms of the shear viscosity, including a shear viscosity which is solely a function 
of temperature, y ~ (mT) 3 / 2 , as predicted by kinetic theory in the dilute limit. 
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I. INTRODUCTION 


Considerable effort has been devoted to extracting transport properties, in particular 
the shear viscosity and the spin diffusion constant, of dilute atomic Fermi gases |l|-|9j. The 
interest in these experiments is driven by the observation that strongly correlated Fermi 
gases can serve as model systems for other quantum many body systems, such as high T c 


superconductors or the quark-gluon plasma 


1CH12||. There are, however, two difficulties that 


have prevented truly model independent measurements of transport coefficients in trapped 
systems so far. The first difficulty is that the diffusion constants for momentum or spin 
depend on the local density while the associated experimental observables are global mea¬ 
sures such as the mean square cloud size or the total spin current. This implies the need to 
unfold the experimental data in order to obtain the density and temperature dependence of 
the transport coefficients. The analogous deconvolution problem for equilibrium quantities 


has been overcome using a number of techniques 


[f! 14], but the first study attempting to 


determine the local shear viscosity only appeared recently [9]. 

The second, more serious, difficulty is that the diffusion approximation breaks down in the 
dilute part of the cloud. This problem cannot be ignored, because a naive application of the 
Navier-Stokes or the diffusion equation to the dilute corona leads to paradoxical behavior. 
Consider, for example, a scale invariant Fermi gas expanding after release from a harmonic 
trap [15]. In the case of a vanishing shear viscosity the expansion dynamics is described by 


an exact scaling solution of the Euler equation. This solution corresponds to a Hubble-like 
flow, in which the fluid velocity u is always linearly proportional to the distance from the 
center of the trap, and the temperature is only a function of time. We can now study how 
this picture is modified in the presence of a small dissipative term. The viscous contribution 
to the stress tensor, 

/ o \ 

( 1 ) 


SH-ij = —r) + VjUi - -SijV ■ u 

is a constant in space that multiplies the local shear viscosity rj. At unitarity scale invariance 
implies that r] = n/(n/T 2 / 3 ), where n is the density, T is the temperature, and f(x ) is a 
universal function. In the dilute limit the shear viscosity is only a function of temperature 
and not of density, 77 = const ■ (mT) 3|/2 . Kinetic theory predicts const = 15/(32i/F) Q, Q. 
Since the temperature is spatially constant one concludes that 5n^ goes to a constant in 
the dilute part of the cloud. This implies that there is no dissipative force, but a constant 
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amount of dissipative heating per unit volume, where the energy for this infinite amount of 
heating is supplied by a heat current that flows in from spatial infinity. 

This description is, of course, completely wrong. The mean free path in the dilute corona 
is much larger than the inter-particle spacing, and there are no collisions that could establish 
dissipative forces or viscous heating. Particles in the dilute corona are ballistically streaming. 
In order to describe the situation correctly, we have to combine a fluid dynamical description 
of the core with a weakly collisional theory of the corona. In this work we suggest that an 
efficient method for achieving this goal is to include certain non-hydrodynamic degrees of 
freedom, an anisotropic pressure tensor, in the theoretical description. In Section [TT] we 
motivate this method by studying certain exact solutions of the Boltzmann equation. In 
Section m we review the derivation of standard, isotropic, fluid dynamics from kinetic 
theory, and in Sections IIVI and [V] we extend this method to anisotropic fluid dynamics. A 
similar method was proposed as an extension of fluid dynamics to describe the early stage 
of a heavy-ion collision, see 18,ll9j. In Sections IVII and [VIII we describe numerical methods 
and show results from an anisotropic fluid dynamics code. This code is a generalization of 
the Navier-Stokes code described in 20 ]. We show that our method describes the crossover 
from fluid dynamics to free streaming both globally, for a shear viscosity of the form r] ~ n, 
and locally, for a shear viscosity of the form 77 ~ (mT) 3//2 . We end with an outlook in 
Section IVIIII 


II. GLOBAL CROSSOVER FROM FLUID DYNAMICS TO BALLISTIC EXPAN¬ 
SION 


The global crossover from fluid dynamics to free streaming in the expansion of a trapped 


Fermi gas after release from the trap was studied in 
to the Boltzmann equation obtained in 22 


2 l|, based on a set of scaling solutions 


23]. We will use these solutions to motivate 


an extension of the fluid dynamic equations that accommodates the transition from fluid 
dynamics to free streaming locally. This a ppr oach is described in Sect. IIVI 

The scaling solutions introduced in [22|, [23l] solve the Boltzmann equation in a harmonic 
confinement potential and using the Bhatnagar-Gross-Krook (BGK) approximation. This 
approximation is based on a collision term of the form C[f] = —Sf/r, where 5f is the 
deviation of the distribution function / from the local equilibrium distribution, and r is the 
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relaxation time. The authors of 2 ll423j further assumed that r is only a function of the local 
temperature, but not of the local density. In the case of two-dimensional traps extensions 
of the scaling ansatz to anharmonic confining potentials were studied in 24], 

Fluid dynamics corresponds to the limit r —> 0, and free streaming is realized as r —> oo. 
In both limits the Boltzmann equation is solved by a distribution function of the form 22] 


/(£, v, t ) = r (t) fo(R(t), U (t)), 


( 2 ) 


where 


MR, U) ~ exp V + Uf] J 

is the initial distribution function in a harmonic potential with frequencies u>i and 


r(i) = n 


„ , . Xi . i Vi — Ui , . bAt) 

Ri(t ) = , Ui{t) = , ai (t) = J 


( 3 ) 


( 4 ) 


^ b^t) m i/2 ’ w bi®’ w 

with Ui = o>i(t)Xi for fixed i. We note that the ansatz in Eq. ([2]) preserves the shape of 
the initial Boltzmann distribution in the cartesian directions i = x,y, z, and that 6{ plays 
the role of an anisotropic scale factor for the temperature. In the free streaming limit the 
solution of the Boltzmann equation is 


0*(f) = , W*) = i 1 + ^t 2 ) 


1/2 


( 5 ) 


In the limit of ideal fluid dynamics we find 6i(t) = 6{t) with 6{t) = [Tl, b t (t)\ 2 / 3 , which 
implies that the temperature is isotropic. The scale parameters bi(t) are determined by 

,2 


bi(t) = 




l 2/3 

n jbj(t) bi(t) 


( 6 ) 


This equation can be solved analytically in the limit of late times and a strongly deformed 
trap, a i z = Acu_l and u> x = u> y = cj± with trap deformation A « 1. In this case one finds 

b±(t) ~ y/3/2 U!j_t. 

Solutions for the transverse flow velocity and the density distribution in the transverse 
plane are shown in Fig. [0 We observe that the free streaming and fluid dynamic solutions 
are qualitatively similar. Transverse pressure in fluid dynamics leads to acceleration, which 
is reflected in the larger expansion velocity of the fluid dynamics solution in the left panel. 
Over time, the larger velocity shifts the peak of the density distribution to larger radii, as 
shown in the right panel. It is interesting to note that the velocity field at late times is the 
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FIG. 1: Comparison between solutions of the Boltzmann equation in the limits of free streaming 
(blue solid curves) and ideal fluid dynamics (red dashed curves). In the left panel we show the 
velocity field component u x and the transverse density profile xn(x) at an early time t = wj 1 . In 
the right panel we show these observables at a later time t = 6WJ 1 . The solutions correspond to 
a trap deformation A = 0.045. The density n and the velocity u x are given in arbitrary units, but 
the scales in the left and the right panel are identical. 

same in free streaming and ideal fluid dynamics. The mean velocity, that is the velocity 
weighted by the density, is larger in fluid dynamics because the maximum of the density 
is shifted to larger radii. In the limit A 1 this difference in the mean velocity can be 
understood in terms of energy conservation. In free streaming the internal energy of the 
fluid is transferred equally to kinetic energy in all three directions. In fluid dynamics most 
of the energy is transferred to transverse motion, and the mean velocity is larger by a factor 
^/2 ~ 1 . 22 . 

Figure [T| shows the difference between ideal fluid dynamics and free streaming in the 
idealized situation that the relaxation time is not a function of density, so that the entire 
cloud is either in the ballistic or the fluid dynamical regime. In reality a transition between 
the two regimes occurs in the dilute part of the cloud, and the transition region may shift 
during the evolution. In Sect. HVI we will study a theoretical approach that can dynamically, 
as a function of time and the spatial coordinates, accommodate the crossover from fluid 
dynamical to ballistic behavior. 
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III. FLUID DYNAMICS FROM KINETIC THEORY 


Before we introduce anisotropic fluid dynamics we review the derivation of standard fluid 
dynamics from kinetic theory. We can view fluid dynamics as an effective description of 
a fluid that arises from the Boltzmann equation in the limit of a short mean free path. 
Consider the Boltzmann equation 




+ v-V x -F-V p )f p (x : t) = C[f p ], 


(7) 


where f p is the single-particle distribution function, C[f p ] is the collision term, v = V p E p the 
velocity of a particle with energy E p , and F = —'V x E p is a force. Using the properties of the 
collision term, in particular the conservation of particle number, energy and momentum, we 
can derive conservation laws for the conserved currents. Taking moments of the Boltzmann 
equation we find (the repeated index j is summed over) 

dp _ 

m +v "* = °- 


dS -* 

IV + V ' J = 0. 


( 8 ) 


The conserved charges, the mass density p = mn, the momentum density If, and the energy 
density £, are given by 

p(x,t) = J dY p rnf p (x, t ), 

If (x,t) = J dT p mvf p (x,t), 

£(x,t) = J dY p E p f p (x, t ), (9) 

where dY p = d 3 p/(2n) 3 . The momentum density n is also the conserved current associated 
with the conservation of mass. The remaining conserved currents are the stress tensor n^- 
and the energy current j £ , 


Ylij(x : t ) = j dT p piVj f p (x, t ), 

f(x,t) = J dVpE p (VpUp) f p (x,t). 


( 10 ) 

( 11 ) 


In order for Eqs. (jHJ) and (J9]) to close we have to supply constitutive equations, that is 
explicit expressions for the conserved currents in terms of the fluid dynamical variables p, n 
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and £. In kinetic theory constitutive equations can be derived by expanding the distribution 
function around the local thermodynamic equilibrium distribution f p , 

/» = f° + Vp + Sf r + ■ ■ ■ , (12) 


where 

/° = exp([p-U(|u-w|)]/T) , (13) 

and 5fp are terms that contain n’th order gradients of the fluid dynamical variables. The 
equilibrium distribution function is expressed in terms of intensive quantities, the local 
chemical potential p(x,t), the temperature T(x,t), and the fluid velocity u(x,t). From 
Eq. (TT3P we can compute the conserved currents at zeroth order in the gradient expansion. 
We get 7 f = pu and 

n ij = puiUj + P5ij , (14) 


as well as j £ = u (w + . Here, P is the pressure and w = £°+P is the enthalpy density, 

where £° denotes the energy density in the local rest frame of the fluid, £° = £ — 4 pu 2 . 
The conservation laws combined with Eq. (fT4]l lead to the Euler equations of ideal fluid 
dynamics. The final ingredient needed to complete the description is an equation of state, 
P = P{£ 0 ,p). Using the dispersion relation of a free particle, E(v) = ^mv 2 , we obtain 
P = -,£()■, which agrees with the exact result for a scale invariant fluid. 

The local equilibrium distribution function is a solution of the Boltzmann equation at 
leading order in the Knudsen number Kn = l m f p /L, where l m f p is the mean free path and L 
is the characteristic distance over which the conserved charges vary. At next order a solution 
can be found most easily by using a very simple form of the collision term. Using the BGK 
collision term 


f - 

c[f p ] = - Jp Jp 


T 


and, again, taking the dispersion relation to be that of a free particle, we find 

'5 T 


Wl_ mT fp 
Jp 2T 


Ci Cj (Jij H - 


L m 


CkQk 


(15) 


(16) 


where we have defined c = v — u, and repeated indices z, j, k are summed over. We have also 
introduced the strain tensor 


(J ij V i Uj 


V j'U'i 




( 17 ) 
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as well as q — — Vlog(T). The corresponding corrections to the conserved currents are 

= —rjCTij , 8jf = UjSUij - acVjT , (18) 

where we have defined the shear viscosity rj = rP and the thermal conductivity k = | rP. 
Incorporating the gradient corrections in Eq. (fT8]l into the conservation laws leads to the 
Navier-Stokes equation. Note that within the approximations used here the shear viscosity 
and the thermal conductivity are proportional to one another, and the bulk viscosity £ is 
zero. In general, rj, ( and k are independent parameters, but in a scale invariant fluid £ = 0 
is an exact result. 


IV. ANISOTROPIC FLUID DYNAMICS FROM KINETIC THEORY 


The gradient expansion fails in the dilute regime of the cloud. An obvious solution to this 
problem is to consider the full Boltzmann equation, see 25]. Such an approach is motivated 
by the observation that even though the classical Boltzmann equation is only justified in the 
dilute regime, it reproduces the fluid dynamical limit in the dense regime. This implies that 
coarse grained observables extracted from the Boltzmann equation are in fact more reliable 
than the kinetic theory which is used to derive them. There are, however, some difficulties 
with this approach. First of all, the Boltzmann equation involves a six-dimensional phase 
space distribution function, and is considerably more difficult to solve than the Navier-Stokes 
equation. Second, the transport properties are now encoded in a non-linear collision integral, 
which is difficult to compute from first principles, and not easy to parameterize in a way 
that allows for a shear viscosity which is a general function of density and temperature. And 
finally, it is difficult to incorporate the empirical equation of state. 

An alternative approach is to use a set of fluid dynamical equations which is equivalent 
to the approach presented in the previous section at some fixed order in the gradient expan¬ 
sion, but also contains extra, non-hydro dynamic, degrees of freedom that ensure a smooth 
crossover to the ballistic regime. Consider 


f P = fr + sf' + Vp + ■ ■ ., 


(19) 


where 


fr = exp 




E 


mc a 

2 T n 


T le = n r 


1/3 


( 20 ) 




The form of the anisotropic distribution function /“ is motivated by the observation that, 
for suitable choices of /i, T a and u a = v a — c a Eq. (1201) is an exact solution of the Boltzmann 
equation describing the expansion from a harmonic trapping potential in the free streaming, 
collisionless limit, see Section EU In order to derive the conservation laws and constitutive 
equations we will use the free dispersion relation E p = p 2 /(2m). This is sufficient in order 
to recover the ballistic and fluid dynamical limits, but restricts the form of the equation 
of state to P = nT = ^Y^ a T a - This is not a problem in the scale invariant limit, because 
the evolution equations are only sensitive to the relation P(£°) = |£°, which is fixed by 
scale invariance. The full equation of state, P = P(n, T), is needed to determine the initial 
density profile from the equation of hydrostatic equilibrium, VP = — nW, where V is the 
confining potential. 

We note that the ansatz for f p n breaks rotational invariance. This particular ansatz is 
intended for analyzing the expansion of a gas cloud from a harmonic confinement potential, 
where the symmetry axes of the potential are aligned with the cartesian coordinate system 
used. Rotational symmetry can be restored by using the more general ansatz 

/„“ = exp [jr-Y. jcJm) , T le = (det (fT 1 ))' 73 . (21) 

For our purposes we will continue to use the simpler ansatz given in Eq. (l20]h 

We can use Eqs. ((9]) - (fTTf to determine the constitutive equations. We find n = pu and 

£ = ipS 2 + f°, £° = \P- (22) 

The stress tensor is given by 

n ij P'Ui'Uj T P ^ij V SUij , <5IRj S m S ja AP a , (23) 

where A P a = P a — P. We use the convention that repeated vector indices i, j, k are summed 
over, but repeated anisotropic indices a, b are not, unless an explicit summation symbol 
occurs. The components of the energy current are 

jf = Ui + wj + Sjf , Sjf = UjSUij = 5 ia u a AP a , (24) 

where w = S° + P. In kinetic theory we also find P = nT and P a = nT a . Combining 
the constitutive equations (T22li - ()24ll with the conservation laws Eqs. (JSJ) and ((9]) gives five 
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equations for seven fluid dynamical variables, /i, P a and iq. We can get two additional 
equations by considering further moments of the Boltzmann equation. The conservation 


laws arise from taking moments with respect to the conserved quantities 1, 

Taking moments with mv 2 /2 (no sum over a) gives 

mv, and mv 2 / 2. 


d£ a - ^ A P a 

3t Ja 2t 

(25) 

where we have defined 

£a — 2^ Ua + ^ a ’ 

(26) 


{fa)i = U i ( 2 ^“ + + ^ia,P^j + (8ja)i > 

(27) 

and £ ° = \P a as well as 

^3a)i ^ia^a^Pa • 

(28) 


Note that Eq. (125]) . when summed over a, gives the equation of energy conservation. The 
remaining two equations determine the non-equilibrium pressure components P a . Also note 
that we can derive additional equations based on the off-diagonal moments with mv a Vb/ 2 
(a 7 ^ b ). These equations determine the off-diagonal components of the temperature in 
Eq. fl2U). 

We will show in Section |V] that Eq. (l23|i reduces to the Navier-Stokes stress tensor in 
the limit r —> 0. This implies that f^ n already contains all terms of order 0(ViUj) in Sfp, 
and that df^ 1 only includes terms associated with heat conduction, 5f = 0(V{T). It is 
straightforward to include these effects, but in the context of an expanding gas cloud heat 
conduction is a very small effect, because the initial state is isothermal and this property 
is preserved by the evolution in ideal fluid dynamics. This implies that gradients of the 
temperature are proportional to r, and, thus, heat flow is a second order effect in the 
relaxation time, acV {T = 0(t 2 ). 

V. FLUID DYNAMICAL EQUATIONS IN LAGRANGIAN FORM 

In practice we solve the equations of fluid dynamics in Lagrangian form. We introduce 
the comoving time derivative Dq = do + u ■ V. The continuity equation can be written as 

Dqp = —pV • u (29) 
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( 30 ) 


and the equations of momentum and energy conservation are 



P 


(31) 


where we have defined the energy per mass e = £/p and 5jf = UjfiUij as v 
5 ia 5j a AP a . The equation for the anisotropic energy density can be written as 


UjdHij as well as <51!^ 



(32) 


where e a = £ a /p and = S ia UjSHij. In standard fluid dynamics we view p, iq and £ 

as the fluid dynamical variables. Their time evolution is governed by Eqs. (J29j) - (j3TT) . and 
in order to determine the RHS of Eqs. (1301) and d3T|) we use the equation of state P(£°) 


with £° = £ — \pu 2 , where in a scale invariant fluid P(£°) = |£°. In anisotropic fluid 


dynamics we have two extra variables, £\ and £2 with = £■ Their time evolution is 

governed by Eq. (l32lh and P a is given by the anisotropic equation of state P a (£°) = 2 £° 
with £® — £„ — \pu 2 a . Note that P = \ J2 a Pa satisfies the isotropic equation of state. 

In the previous section we argued that in the limit r —)• 0 anisotropic fluid dynamics re¬ 
duces to Navier-Stokes viscous fluid dynamics. We expect, in particular, that the dissipative 
correction to the stress tensor hfl^ = 5i a Sj a AP a approaches 8Ti l3 = —r)o %3 (for i = j) with 
p = tP. To see this we rewrite Eq. (1321) as 



(33) 


and solve for AP a at leading order in r. This implies that in evaluating £® and (<5jf)* we 
can replace P a by P, so that e a = — |n 2 + \u 2 a and (5jf), = 0. We use the equations of 

ideal fluid dynamics to compute D 0 e and D 0 Ui and find 



(34) 


This result shows that anisotropic fluid dynamics relaxes to the Navier-Stokes equation with 


(35) 


V = tP , 


where r is the relaxation time. Note that the expression in Eq. (f23ll does not reproduce the 
off-diagonal components of in Navier-Stokes theory. In order to study flows in which 


these terms are non-zero we have to start from the more general ansatz in Eq. (12T| and 
consider moments of the Boltzmann equation with mv a Vb/2 for a 7 ^ b. 
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VI. ANISOTROPIC FLUID DYNAMICS: NUMERICAL METHOD AND 
CHOICE OF UNITS 


We have implemented anisotropic fluid dynamics as an extension of the Navier-Stokes 


code described in 


20]. The Navier-Stokes code solves the advection equations in Lagrangian 


coordinates. A Lagrangian time step is followed by a piecewise parabolic remap onto 
an Eulerian grid. The algorithm is based on the PPMLR (Piecewise-Parabolic Method, 
Lagrangian-Remap) scheme developed by Colella and Woodward |26j and implemented as a 
multi-dimensional method in the VH1 code written by Blondin and Lufkin 27]. The main 
modification is that we add the fluid dynamical variables P a and £ a and solve Eq. (TT?1) . We 


solve for all three components of P a and verify that | ]T a Pa agrees with P. Since Eq. 

is a relaxation equation for A P a we have to choose the time step as 

. . ( Ax Ax 

At = mm c s -, Ui -, r 

x V 2 2 ’ . 


(36) 


where c s is the local speed of sound, Ax is the grid spacing, and r = r\/P is the local relax¬ 
ation time. The first two constraints arise from the condition that disturbances emerging 
from opposite faces of a fluid cell cannot interact during a time step. The third constraint 
ensures that the relaxation time equation is stable. The condition At <r implies that the 
simulation becomes inefficient if r is very small, which is close to the limit of ideal fluid 
dynamics. In principle this can be addressed by using the analytic result given in Eq. (j34lh 
but we have not done so in the present work. 

We have studied the evolution of a unitary Fermi gas after release from a harmonic trap. 
The trapping potential is V(x ) = ^mufx f with u x = u y = u± and cm = Acnj_- We use 
dimensionless variables for distance, time and velocity based on the following system of 
units 20] 

p/6 ( 


xq = (3N\y 


\3mu 




1/2 


to — aR 1 ) 


u o = x 0 u± 


(37) 


The unit of density is n 0 = x 0 3 . The corresponding units for energy density, pressure, and 
temperature are given by 


£o = 


mu;, 


Xq 


Pq = 


mu i 


Xq 


Tq = mu , x ; 


2 J2. 


jao • 


Finally, the unit of the shear viscosity is 


Vo = 


mu± 

Xq 


(38) 


(39) 
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In the high temperature limit the initial density is a Gaussian. The central density is given 
by 

where N is the number of particles, Ep = (SiVA) 1 / 3 ^ is the Fermi energy, and Eq is the 
total energy per particle of the trapped gas. Moreover, it is convenient to normalize the 
dimensionless central density n(0)/n 0 to one. This means we also divide the density by 
the dimensionless factor ( NX)/tt 3 / 2 • (Ep/E 0 ) 3 / 2 in equ. (1401) . The normalized dimensionless 
shear viscosity is 

__ V n o 
V Von(0)' 

In the following we will consider a shear viscosity of the form r] = a n n + ap(mT) 3 ^ 2 . The 
corresponding dimensionless shear viscosity is fj = a n n + apT 3 ^ 2 with 

4tt 3 / 2 a T (E 0 \ 3/2 


Oi r , 


Otr). 


2 (31VA) 1 / 3 ’ 3 (3N\y/ 3 \EpJ 

Kinetic theory predicts that in the high temperature limit a n = 0 and ap = 15/(32y / 7r). In 
the anisotropic fluid dynamics framework the shear viscosity is determined by the relaxation 
time. The dimensionless relaxation time is 


(42) 


r r) 

T = h = p- 


(43) 


Equation (I42]l shows that, in dimensionless units, the number of particles N only appears 
together with the viscosity coefficient. This implies that the ideal evolution is independent 
of N, and that for a given viscosity dissipative effects are smaller for a larger number of 
particles. 


VII. ANISOTROPIC FLUID DYNAMICS: RESULTS 

We first consider the case ap = 0 and study the dependence on a n . Figure [2] shows 
the time evolution of the aspect ratio Ap(t) = (rp)/(r 2 ), defined by the ratio of mean 
squared transverse and longitudinal cloud radii, for different values of the shear viscosity, 
a n = 0.1, 1, 1000. For comparison we also show the result in ideal fluid dynamics, the free 
streaming limit, and the solution of the Navier-Stokes equation for a n = 0.1, 1. We consider 
a Gaussian initial condition, which corresponds to a solution of the hydrostatic equation in 
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FIG. 2: This figure shows the evolution of the aspect ratio Ar as a function of time t in dimension¬ 
less units, as explained in the text. The initial trap deformation is Afl(O) = A = 0.1 and the initial 
energy is Eq/E p = 1. The red solid curve shows the evolution in ideal fluid dynamics, and the 
black dashed curve is the free streaming limit. The remaining curves were obtained using viscous 
fluid dynamics with different values of a n for the shear viscosity fj = a n fi. The dotted curves show 
Navier-Stokes results for a n = 0.1, 1 and the data points are the corresponding predictions from 
anisotropic fluid dynamics. In the case of anisotropic fluid dynamics we also show the result for 
a n = 1000, which is close to the free streaming limit. 

the case of an equation of state of a free gas, P = nT. At t = 0 the aspect ratio is given by 
the trap deformation, A/j(0) = A. Pressure gradients preferentially accelerate the fluid in the 
transverse direction, and An(t) grows as a function of time. Viscous effects counteract the 
expansion in the transverse direction, and accelerate the fluid in the longitudinal direction, 
reducing the value of An{t). 

For the smallest value of the shear viscosity, a n = 0.1, we find good agreement between 
anisotropic fluid dynamics and Navier-Stokes theory, as expected from Eq. (l34]h For larger 
values of a n anisotropic fluid dynamics predicts that dissipative effects saturate, and that 
A R (t) approaches the free streaming limit. In Navier-Stokes theory, on the other hand, 
dissipative effects continue to grow with a n and the evolution of A R (t) becomes arbitrarily 
slow. 

More details are provided by Fig. [3l In this figure, we compare the dissipative corrections 
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x/x 0 


FIG. 3: This figure shows the xx component of the dissipative correction to the stress tensor 
8H xx (x, 0,0) in Navier-Stokes theory (green solid curve) and in anisotropic fluid dynamics (red 
squares) for different times t/to = 0.75 — 2.0 in steps of At = 0.25to- The magnitude of the 
viscous stresses decreases with time. The calculation was performed with a density dependent 
shear viscosity fj = a n fi and a n = 0.15. We used a trap deformation A = 0.045 and an initial 
energy E 0 /E F = 3. 


to the stress tensor in Navier-Stokes theory and in anisotropic fluid dynamics. We focus 
on the xx components = —rja xx and 5U XX = A P x at different times during the 

evolution of the expanding gas cloud. The calculation was performed for a n = 0.15, so that 
the aspect ratio A F (t) shows good agreement between Navier-Stokes theory and anisotropic 
fluid dynamics. Note that SU^x hi Fig. [3]was computed using the velocity field in anisotropic 
fluid dynamics. We observe that the two dissipative stress tensors are indeed very close, and 
that the agreement improves at late times. This indicates that the equations of anisotropic 
fluid dynamics contain second order terms in r that describe the relaxation of the stress 
tensor to the Navier-Stokes limit 
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29]. 


Figure [4] demonstrates that anisotropic fluid dynamics can be applied to the case of a 
purely temperature dependent shear viscosity, r] = « T (mT) 3 ^ 2 , for which Navier-Stokes fluid 
dynamics fails. We observe that in the center of the cloud the two dissipative corrections 
to the stress tensor are close, in particular at late times. In the corona, however, SH^ 


and are very different. As explained in Section [D the dissipative contribution to the 
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FIG. 4: Same as Fig. [3] for a temperature dependent shear viscosity 77 = a^F 3 / 2 with ay = 0.06. 
The blue dashed line shows the negative of the anisotropic pressure component P x . 


Navier-Stokes stress tensor is approximately constant in space. In contrast, the dissipative 
contribution to the stress tensor in anisotropic fluid dynamics goes to zero in the dilute part 
of the cloud. As we can see from the blue dashed curves in Fig. 0] this happens in the regime 
where the dissipative stresses are comparable to the total pressure of the fluid, |5II^f | ~ P x . 
This condition, corresponding to the point where the dashed blue line intersects the red 
symbols, signals the breakdown of Navier-Stokes theory. 

We note that — W X 8U XX corresponds to a force that points towards the center of the 
cloud, and reduces the transverse expansion. We also observe that at late times a xx , which 
measures the slope of the velocity field, is smaller in the center than in the corona. Viscous 
forces push on the center of the cloud and slow it down, while the ballistic corona is lifting 
off. 


In anisotropic fluid dynamics the viscous stresses are concentrated in the center of the 
gas cloud even if rja i3 is not localized. We therefore expect that the time evolution of the 
aspect ratio A R (t) can be described by an effective density dependent shear viscosity 77 ~ n 
even if the microscopic shear viscosity is only a function of temperature. This approximation 
has been used to analyze experimental data on the expansion of trapped Fermi gases near 
unitarity jd, 8 ]. In Fig. [5] we show that for a given initial temperature, and for a suitably 
chosen value of a n , the evolution of A R (t ) is indeed essentially indistinguishable between 
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FIG. 5: This figure shows the time evolution of the aspect ratio Ar from a simulation using 
anisotropic fluid dynamics with = 0.023 in 77 = d^T 3 / 2 and A = 0.1, Eq/Er = 1 (blue squares). 
For comparison we also show the result in ideal fluid dynamics (red solid curve) and the case 
of purely ballistic expansion (black dashed curve). The green dotted curve shows a fit to the 
ax = 0.023 result based on Navier-Stokes theory with a n = 0.066 in fj = a n n. 


the two cases r/ ~ n and 77 ~ [mTf/ 2 . 

In order to resolve this ambiguity and determine the full microscopic dependence of 77 
on n and T we have to study the sensitivity of the effective a n on the initial temperature 
T or the initial energy Eq/Er. In the literature, the value of a n which describes the time 
evolution of Aji(t) at a given initial temperature is referred to as the trap averaged value 
of 77/77., denoted (a n ) |4|. In Fig. [ 6 ] we show the dependence of (a n ) on E 0 /E F for the 
case 77 = 15/(32^/77) (777T) 3 / 2 . We consider a Gaussian initial condition, so that E 0 = 3 T. 
We observe that the growth of ( a n ) is not simply proportional to E^ 2 . This is because 
the evolution is not only sensitive to the temperature dependence of the shear viscosity, 
77 ~ T 3 / 2 , but also to the temperature dependence of the relaxation time, r ~ T 1 / 2 , and the 
temperature dependence of the effective relaxation volume. 

We note that as a consequence of the complicated dependence of ( a n ) on the system size 
and lifetime the result is not universal, which means that ( a n ) depends on the number of 
particles N and the trap deformation A. In the present work we will not attempt to perform a 
detailed analysis of the experimental data obtained in j^, 9], This will require implementing 
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FIG. 6: This figure shows the dependence of (a n ), the effective trap averaged ratio 77 /n, on the 
initial energy per particle in units of the Fermi energy Ep. A precise definition of (a n ) is given 
in the text. The calculation was performed for a gas cloud of 2 • 10 5 particles with an initial trap 
deformation of A = 0.045. The microscopic shear viscosity is given by the kinetic theory result in 
the high temperature limit, 77 = 15/(32- v /7r)(mT) 3 / 2 . 


a non-axially symmetric confining potential, a non-Gaussian initial density distribution, a 
realistic equation of state P(n,T ), and a shear viscosity which is a function of both n and 
T. It is nevertheless interesting to consider a rough comparison between our results and the 
analysis in [9]. At E 0 /E F = 3.14 we find ( a n ) = 13.49, compared to (a n ) = 19.63 ± 0.54 
reported in [9]. The discrepancy is an indication of the magnitude of the effects due to the 
trap geometry, contributions beyond the dilute limit of the equation of state and the shear 
viscosity, and possible shortcomings of our description in the regime where the size of the 
fluid dynamical core shrinks to zero, and a kinetic treatment of the entire cloud is more 
appropriate. 


VIII. CONCLUSIONS AND OUTLOOK 

In this work we have shown that by including non-hydrodynamic degrees of freedom it 
is possible to achieve a smooth transition between fluid dynamics and the ballistic limit. 
There are a number of interesting applications and theoretical issues that remain to be 
investigated. 
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• Other fluid dynamic problems: The anisotropic fluid dynamics approach, described by 
equ. (I291l32p . provides a general scheme for addressing problems in which the transition 
regime from fluid dynamics to ballistic behavior plays a role. In this work we have 
applied this method to the expansion of a unitary Fermi gas from a harmonic trap, but 
the applicability of the approach is clearly much broader, including Bose or classical 
gases, as well as other experimental observables, such as collective modes. 

• Restoring rotational invariance: In this work we only considered the “cartesian” ansatz 
in Eq. (l20j) . In order to restore full rotational invariance we have to start from Eq. (12T1) 
and derive the corresponding fluid dynamical equations. 


• Relation to second order fluid dynamics: In our numerical simulations we observed 
that anisotropic fluid dynamics contains some effects that appear at second order in 
the gradient expansion of fluid dynamics, in particular a finite viscous relaxation time. 
It will be interesting to make this more precise, and extend the equations of motion 
to complete second order accuracy. 


More accurate treatment of the Knudsen limit: We have shown that anisotropic fluid 
dynamics reproduces the Navier-Stokes equation at order 0(r ), where r is the relax¬ 
ation time in the BGK approximation to the Boltzmann equation. In the opposite 
limit, r —> oo, anisotropic fluid dynamics provides an exact solution of the Boltzmann 
equation for r = oo. However, at 0{l/r) anisotropic fluid dynamics does not provide 
an exact solution of the Boltzmann equation, only a solution at second order in an 
expansion in moments of the distribution function with respect to momentum. The 
reliability of this approximation can be studied by comparing with the exact numerical 
solutions obtained in 25]. 

Extension to superfluid hydro dynamics: The experimental results obtained in 9] also 
cover the regime T < T C1 where T c is the critical temperature for superfluidity. In 
order to extract the shear viscosity in this regime we have to extend anisotropic fluid 
dynamics to the superfluid (two-fluid) regime. 


• Extension to spin diffusion: The problem related to the dilute regime also affects the 
extraction of the spin diffusion constant. It will be interesting to study whether our 
method can be extended to the case of charge and spin diffusion. 
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Our immediate goal is to use the method developed in this work to extract the local 
shear viscosity rj(n,T) = nf{n/T 2 / 3 ) from the data presented in |q]. This can be achieved 
by inverting the dependence of the aspect ratio A R (E 0 ,t) as a function of the initial energy 
on the function f(x). As explained in the previous section, this will require implementing 
a non-Gaussian initial density distribution as well as considering a non-axially symmetric 
trap, and a realistic equation of state P(n, T ). A natural starting point for unfolding the full 
density and temperature dependence of the shear viscosity is the reconstruction presented 
in { 9 ]. At this point we have not implemented a superfluid version of the anisotropic fluid 
dynamics method, and we are limited to the regime T > T c . 
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